##########################################################################################

library(data.table)
library(optparse)
library(ggplotify)
library(ggsci)

##########################################################################################
option_list <- list(
    make_option(c("--input_file"), type = "character"),
    make_option(c("--scriptPath"), type = "character"),
    make_option(c("--out_path"), type = "character") 
)

if(1!=1){
    
    ## 活性peak和motif的数量文件
    input_file <- "~/20231121_singleMuti/results/motif_position/All.motif_peak.tsv"

    ## 既往研究整理的代码
    scriptPath <- "~/20231121_singleMuti/scripts/scScalpChromatin"

    ## 输出
    out_path <- "~/20231121_singleMuti/results/motif_position"

}

###########################################################################################
parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

input_file <- opt$input_file
scriptPath <- opt$scriptPath
out_path <- opt$out_path

dir.create(out_path , recursive = T)

###########################################################################################
## 导入数据
data <- fread(input_file)

source(paste0(scriptPath, "/plotting_config.R"))

###########################################################################################
## 细胞类型
dodge_width <- 0.75
dodge <- position_dodge(width=dodge_width)

#细胞上色
use_colors <- c(pal_npg("nrc")(10) , pal_jco("default")(6))
#names(use_colors) <- unique(scrnat$cell_type)
names(use_colors) <-c( "Myoid cells","Leydig cells" ,          
"Endothelial cells","Zygotene",         
"Round&ElongateS.tids","Patchytene",
"SSC","Sperm" ,                 
"Diplotene","Early stage of spermatids", 
"Leptotene","Sertoli cells",            
"Macrophages","Differenting&Differented SPG",
"Pericytes","NKT cells")

## 细胞顺序
grp_order = c("SSC", "Differenting&Differented SPG", "Leptotene",
"Zygotene", "Patchytene", "Diplotene",
"Early stage of spermatids", "Round&ElongateS.tids", "Sperm",
"Leydig cells", "Myoid cells", "Pericytes",
"Sertoli cells", "Endothelial cells", "NKT cells", "Macrophages")

###########################################################################################
## peak
#peak_num（log10转化）

width <- 10
height <- 9

data$cell_type <- factor(data$cell_type, levels = grp_order)
data$log10peak_num <- log10(data$peak_num)
p1<- (
    ggplot(data, aes(x=cell_type, y=log10peak_num, fill=cell_type))
    + geom_violin(aes(fill=cell_type), alpha=0.5, adjust = 1.0, scale='width', position=dodge)
    + geom_boxplot(alpha=0.5, width=0.25, outlier.shape = NA)
    + scale_color_manual(values=use_colors, limits=names(use_colors), name="cell_type", na.value="grey")
    + scale_fill_manual(values=use_colors)
    + guides(fill=guide_legend(title="cell_type"), 
      colour=guide_legend(override.aes = list(size=5)))
    + xlab("")
    + ylab("Peak num(log10)")
    + theme_BOR(border=FALSE)
    + theme(panel.grid.major=element_blank(), 
            panel.grid.minor= element_blank(), 
            plot.margin = unit(c(0.25,1,0.25,1), "cm"), 
            #aspect.ratio = aspectRatio, # What is the best aspect ratio for this chart?
            legend.position = "none", # Remove legend
            axis.text.x = element_text(angle = 90, hjust = 1))
    + scale_y_continuous(limits=c(2,5), expand = c(0, 0))
)

out_file <- paste0(out_path, "/scATAC_peak_num_per_cell_violin1.pdf")
ggsave(  out_file , p1,width=width, height=height)

#motif_num

p2<- (
    ggplot(data, aes(x=cell_type, y=motif_num, fill=cell_type))
    + geom_violin(aes(fill=cell_type), alpha=0.5, adjust = 1.0, scale='width', position=dodge)
    + geom_boxplot(alpha=0.5, width=0.25, outlier.shape = NA)
    + scale_color_manual(values=use_colors, limits=names(use_colors), name="cell_type", na.value="grey")
    + scale_fill_manual(values=use_colors)
    + guides(fill=guide_legend(title="cell_type"), 
      colour=guide_legend(override.aes = list(size=5)))
    + xlab("")
    + ylab("Motif num")
    + theme_BOR(border=FALSE)
    + theme(panel.grid.major=element_blank(), 
            panel.grid.minor= element_blank(), 
            plot.margin = unit(c(0.25,1,0.25,1), "cm"), 
            #aspect.ratio = aspectRatio, # What is the best aspect ratio for this chart?
            legend.position = "none", # Remove legend
            axis.text.x = element_text(angle = 90, hjust = 1))
    + scale_y_continuous(limits=c(200,750), expand = c(0, 0))
)

out_file <- paste0(out_path, "/scATAC_motif_num_per_cell_violin1.pdf")
ggsave(  out_file , p2,width=width, height=height)
